#include "fluid.h"
#include <fstream>
#include <sstream>
#include <stdexcept>
#include <algorithm>
#include <cmath>
#include <iostream>

Fluid::Fluid(const std::string& pvtFile) {
    loadPVTFile(pvtFile);

    // 初始化临界参数（以天然气为例）
    criticalPressure_ = 4.6e6;     // 4.6 MPa
    criticalTemperature_ = 190.6;  // 190.6 K
}

void Fluid::loadPVTFile(const std::string& filename) {
    std::ifstream file(filename);
    if (!file.is_open()) {
        throw std::runtime_error("Cannot open PVT file: " + filename);
    }

    std::string line;
    bool readingComponents = false;
    bool readingPVTPoints = false;

    while (std::getline(file, line)) {
        // 跳过空行
        if (line.empty()) continue;

        if (line.find("COMPONENTS = ") != std::string::npos) {
            readingComponents = true;
            continue;
        }

        if (line.find("PVTTABLE POINT = ") != std::string::npos) {
            readingPVTPoints = true;
            
            // 解析PVT点数据
            size_t start = line.find("(") + 1;
            size_t end = line.find(")");
            std::string data = line.substr(start, end - start);
            
            std::stringstream ss(data);
            std::string value;
            PVTPoint point;
            
            // 读取压力、温度和密度数据
            int valueIndex = 0;
            while (std::getline(ss, value, ',')) {
                double val = std::stod(value);
                switch(valueIndex) {
                    case 0: point.pressure = val; break;      // Pa
                    case 1: point.temperature = val + 273.15; // 转换为K
                    case 2: point.gasDensity = val; break;    // kg/m3
                    case 3: point.oilDensity = val; break;    // kg/m3
                    case 4: point.waterDensity = val; break;  // kg/m3
                }
                valueIndex++;
            }
            
            // 添加数据有效性检查
            if (point.pressure <= 0 || point.temperature <= 0 || 
                point.gasDensity <= 0 || point.oilDensity <= 0) {
                std::cout << "Warning: Invalid PVT data point: P=" << point.pressure 
                         << ", T=" << point.temperature << ", rho_g=" << point.gasDensity 
                         << ", rho_o=" << point.oilDensity << std::endl;
                continue;  // 跳过无效数据点
            }
            
            pvtData_.push_back(point);
        }
    }

    // 打印读取的数据点数量
    std::cout << "Loaded " << pvtData_.size() << " PVT data points" << std::endl;
    if (!pvtData_.empty()) {
        std::cout << "P range: " << pvtData_.front().pressure << " - " 
                 << pvtData_.back().pressure << " Pa" << std::endl;
        std::cout << "T range: " << pvtData_.front().temperature << " - " 
                 << pvtData_.back().temperature << " K" << std::endl;
    }

    if (pvtData_.empty()) {
        throw std::runtime_error("No PVT data points found in file");
    }
}

double Fluid::interpolate(double p, double T, int propertyIndex) const {
    // 检查输入值是否在有效范围内
    if (pvtData_.empty()) {
        return 0.0;  // 如果没有PVT数据,返回0
    }

    // 找到压力和温度的范围
    double minP = pvtData_[0].pressure;
    double maxP = pvtData_.back().pressure;
    double minT = pvtData_[0].temperature;
    double maxT = pvtData_.back().temperature;

    // 如果超出范围,使用边界值
    p = std::max(minP, std::min(p, maxP));
    T = std::max(minT, std::min(T, maxT));

    // 找到压力和温度的插值区间
    auto lowerP = std::lower_bound(pvtData_.begin(), pvtData_.end(), p,
        [](const PVTPoint& point, double val) { return point.pressure < val; });
    auto lowerT = std::lower_bound(pvtData_.begin(), pvtData_.end(), T,
        [](const PVTPoint& point, double val) { return point.temperature < val; });

    // 如果找不到合适的区间,使用最近的值
    if (lowerP == pvtData_.begin()) lowerP++;
    if (lowerT == pvtData_.begin()) lowerT++;
    if (lowerP == pvtData_.end()) lowerP--;
    if (lowerT == pvtData_.end()) lowerT--;

    const PVTPoint& p1 = *(lowerP - 1);
    const PVTPoint& p2 = *lowerP;
    const PVTPoint& t1 = *(lowerT - 1);
    const PVTPoint& t2 = *lowerT;

    // 计算插值系数
    double fp = (p - p1.pressure) / (p2.pressure - p1.pressure);
    double ft = (T - t1.temperature) / (t2.temperature - t1.temperature);

    // 根据属性索引选择要插值的物性
    double value = 0.0;
    switch(propertyIndex) {
        case 0: // 气相密度
            value = (1-fp)*(1-ft)*p1.gasDensity + fp*(1-ft)*p2.gasDensity +
                   (1-fp)*ft*t1.gasDensity + fp*ft*t2.gasDensity;
            break;
        case 1: // 油相密度
            value = (1-fp)*(1-ft)*p1.oilDensity + fp*(1-ft)*p2.oilDensity +
                   (1-fp)*ft*t1.oilDensity + fp*ft*t2.oilDensity;
            break;
        case 2: // 水相密度
            value = (1-fp)*(1-ft)*p1.waterDensity + fp*(1-ft)*p2.waterDensity +
                   (1-fp)*ft*t1.waterDensity + fp*ft*t2.waterDensity;
            break;
    }

    // 确保返回值是有效的
    if (!std::isfinite(value)) {
        // 如果插值结果无效,返回最近的有效值
        value = p1.gasDensity;  // 使用参考点的值
    }

    return value;
}

double Fluid::getGasDensity(double p, double T) const {
    return interpolate(p, T, 0);
}

double Fluid::getOilDensity(double p, double T) const {
    return interpolate(p, T, 1);
}

double Fluid::getWaterDensity(double p, double T) const {
    return interpolate(p, T, 2);
}

// TODO: 实现粘度计算函数
double Fluid::getGasViscosity(double p, double T) const {
    // 临时返回一个常数值，后续需要根据实际数据计算
    return 1.8e-5;  // Pa·s
}

double Fluid::getOilViscosity(double p, double T) const {
    // 临时返回一个常数值，后续需要根据实际数据计算
    return 5.0e-3;  // Pa·s
}

double Fluid::getWaterViscosity(double p, double T) const {
    // 临时返回一个常数值，后续需要根据实际数据计算
    return 1.0e-3;  // Pa·s
}

double Fluid::getSpecificHeat(double p, double T) const {
    // 简化模型：假设比热容是常数
    return 2000.0;  // J/(kg·K)
}

double Fluid::getLatentHeat(double p, double T) const {
    // 简化模型：假设汽化潜热是常数
    return 2.0e6;   // J/kg
}

double Fluid::getGasEnthalpy(double p, double T) const {
    // 简化模型：假设理想气体
    double cp = getGasHeatCapacity(p, T);
    return cp * (T - 273.15);  // 参考温度取273.15K
}

double Fluid::getLiquidEnthalpy(double p, double T) const {
    // 简化模型：假设不可压缩液体
    double cp = getLiquidHeatCapacity(p, T);
    return cp * (T - 273.15);
}

double Fluid::getGasHeatCapacity(double p, double T) const {
    // 简化模型：假设常数比热容
    return 2000.0;  // J/(kg·K)
}

double Fluid::getLiquidHeatCapacity(double p, double T) const {
    // 简化模型：假设常数比热容
    return 4000.0;  // J/(kg·K)
} 